rm(list = ls())
require(stats4)

#multinomial coefficient

multinom=function(vec){
  if(length(vec)==2){
    return(sum(log(1:max(1,sum(vec))))-sum(log(c(1:max(1,vec[1]),1:max(1,vec[2])))))
  }
  if(length(vec)==3){
    return(sum(log(1:max(1,sum(vec))))-sum(log(c(1:max(1,vec[1]),1:max(1,vec[2]),1:max(1,vec[3])))))
  }
}

#Data
#Load in some Data
data2=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_may11_type2clean.csv", header=TRUE)
data3=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_may11_type3clean.csv", header=TRUE)
data4=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_may11_type4clean.csv", header=TRUE)
data5=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_may11_type5clean.csv", header=TRUE)

data16=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Columbia_experiment_September27_type16clean.csv", header=TRUE)
data17=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Columbia_experiment_September27_type17clean.csv", header=TRUE)
data18=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Columbia_experiment_September27_type18clean.csv", header=TRUE)
data19=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/Columbia_experiment_September27_type19clean.csv", header=TRUE)

dataframe16=data.frame(data16)
dataframe16=subset(dataframe16,uid>1600)
dataframe17=data.frame(data17)
dataframe17=subset(dataframe17,uid>1700)
dataframe18=data.frame(data18)
dataframe18=subset(dataframe18,uid>1800)
dataframe19=data.frame(data19)
dataframe19=subset(dataframe19,uid>1900)

dataframe2=rbind(data.frame(data2),dataframe16)
dataframe3=rbind(data.frame(data3),dataframe17)
dataframe4=rbind(data.frame(data4),dataframe18)
dataframe5=rbind(data.frame(data5),dataframe19)

#Properly assign parameters to different experiment types
dataframe2$version=1
dataframe2$Bcorrinc=55
dataframe2$Bfailinc=40

dataframe3$version=2
dataframe3$Bcorrinc=52
dataframe3$Bfailinc=40

dataframe4$version=3
dataframe4$Bcorrinc=52
dataframe4$Bfailinc=30

dataframe5$version=4
dataframe5$Bcorrinc=55
dataframe5$Bfailinc=30

#Combine data and assign useful dummies
dataframe=rbind(dataframe2,dataframe3,dataframe4,dataframe5)
dataframe$Bresponse=(dataframe$response==19)
dataframe$Bstate=(dataframe$state==4)
dataframe$Cpresent=(dataframe$qtype==4)

#Show N for each experiment type
#**
length(unique(dataframe2$uid))
length(unique(dataframe3$uid))
length(unique(dataframe4$uid))
length(unique(dataframe5$uid))
#**


#parameters for different treatments
b1vec=c(40,40,30,30)
b2vec=c(55,52,55,52)

#reward level (note that payoffs are in probability points)
r=0.4

#Initialize relevant data moments with C
data_a1_wC=vector("numeric")
data_a2_wC=vector("numeric")
data_b1_wC=vector("numeric")
data_b2_wC=vector("numeric")
data_c1_wC=vector("numeric")
data_c2_wC=vector("numeric")

#Initialize relevant data moments no C
data_a1_noC=vector("numeric")
data_a2_noC=vector("numeric")
data_b1_noC=vector("numeric")
data_b2_noC=vector("numeric")



uidlist=unique(dataframe$uid)
indioutframe=data.frame(uid=vector('numeric'),likelin=vector('numeric'),errlin=vector('numeric'),likepow=vector('numeric'),errpow=vector('numeric'),likegen=vector('numeric'),errgen=vector('numeric'),likeinat=vector('numeric'),errinat=vector('numeric'))
for(indii in 1:length(uidlist)){
indiframe=subset(dataframe, uid==uidlist[indii])


#fill data values
for(expi in 1:4){
playframe=subset(indiframe,version==expi)
playframe_noC=subset(playframe,Cpresent==FALSE)
data_a1_noC[expi]=sum((playframe_noC$response==18)*(playframe_noC$state==3))
data_a2_noC[expi]=sum((playframe_noC$response==18)*(playframe_noC$state==4))
data_b1_noC[expi]=sum((playframe_noC$response==19)*(playframe_noC$state==3))
data_b2_noC[expi]=sum((playframe_noC$response==19)*(playframe_noC$state==4))

playframe_wC=subset(playframe,Cpresent==TRUE)
data_a1_wC[expi]=sum((playframe_wC$response==18)*(playframe_wC$state==3))
data_a2_wC[expi]=sum((playframe_wC$response==18)*(playframe_wC$state==4))
data_b1_wC[expi]=sum((playframe_wC$response==19)*(playframe_wC$state==3))
data_b2_wC[expi]=sum((playframe_wC$response==19)*(playframe_wC$state==4))
data_c1_wC[expi]=sum((playframe_wC$response==22)*(playframe_wC$state==3))
data_c2_wC[expi]=sum((playframe_wC$response==22)*(playframe_wC$state==4))
}

conmat_wC=matrix(c(1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1, -1,-1,0,0, 0,0,-1,-1), nrow=6, byrow=TRUE)
convec_wC=c(0.001,0.001,0.001,0,-0.99,-0.99 )

conmat_noC=matrix(c(1,0, 0,1,  -1,0, 0,-1), nrow=4, byrow=TRUE)
convec_noC=c(0.001,0.001,-0.99,-0.99 )

const_wCvec=vector("numeric")
const_noCvec=vector("numeric")
for(expi in 1:4){
const_wCvec[expi]=(multinom(c(data_a2_wC[expi],data_b2_wC[expi],data_c2_wC[expi])))+(multinom(c(data_a1_wC[expi],data_b1_wC[expi],data_c1_wC[expi])))
const_noCvec[expi]=(multinom(c(data_a2_noC[expi],data_b2_noC[expi])))+(multinom(c(data_a1_noC[expi],data_b1_noC[expi])))
}

const_wC=sum(const_wCvec)
const_noC=sum(const_noCvec)


#################################################
#Linear Mutual Information

choiceprobslin=function(kappa,treatment){
b1pay=b1vec[treatment]
b2pay=b2vec[treatment]
#Without C
#Information cost function
#cprobs in format a1,a2
cost_noC=function(cprobs){
overallcost=kappa/2*(cprobs[1]*log(cprobs[1]/(cprobs[1]+cprobs[2]))+(1-cprobs[1])*log((1-cprobs[1])/(2-cprobs[1]-cprobs[2]))+cprobs[2]*log(cprobs[2]/(cprobs[1]+cprobs[2]))+(1-cprobs[2])*log((1-cprobs[2])/(2-cprobs[1]-cprobs[2])))
return(overallcost)
}

#Optimization target
negpayoff_noC=function(cprobs){
return(-1*(r*(50/2*(cprobs[1]+cprobs[2])+b1pay/2*(1-cprobs[1])+b2pay/2*(1-cprobs[2]))-cost_noC(cprobs)))
}

#Find optimal accuracies
optimout_noC=constrOptim(rep(0.5,2),negpayoff_noC, grad=NULL,ui=conmat_noC,ci=convec_noC)
cprobvec_noC=optimout_noC$par

#With C
#Information cost function
#cprobs in format a1,b1,a2,b2
cost_wC=function(cprobs){
overallcost=kappa/2*(cprobs[1]*log(cprobs[1]/(cprobs[1]+cprobs[3]))+cprobs[2]*log(cprobs[2]/(cprobs[2]+cprobs[4]))+(1-cprobs[1]-cprobs[2])*log((1-cprobs[1]-cprobs[2])/(2-sum(cprobs)))+cprobs[3]*log(cprobs[3]/(cprobs[1]+cprobs[3]))+cprobs[4]*log(cprobs[4]/(cprobs[2]+cprobs[4]))+(1-cprobs[3]-cprobs[4])*log((1-cprobs[3]-cprobs[4])/(2-sum(cprobs))))
return(overallcost)
}

#Optimization target
negpayoff_wC=function(cprobs){
return(-1*(r*(50/2*(cprobs[1]+cprobs[3])+100/2*(1-cprobs[1]-cprobs[2])+b1pay/2*(cprobs[2])+b2pay/2*(cprobs[4]))-cost_wC(cprobs)))
}

#Find optimal accuracies
optimout_wC=constrOptim(rep(0.3,4),negpayoff_wC, grad=NULL,ui=conmat_wC, ci=convec_wC,)
cprobvec_wC=optimout_wC$par

return(list(cprobvec_noC,cprobvec_wC))

}

likelihoodlin=function(kappa,err){
cprobvec_noC=matrix(c(choiceprobslin(kappa,1)[[1]],choiceprobslin(kappa,2)[[1]],choiceprobslin(kappa,3)[[1]],choiceprobslin(kappa,4)[[1]]),ncol=4)
cprobvec_wC=matrix(c(choiceprobslin(kappa,1)[[2]],choiceprobslin(kappa,2)[[2]],choiceprobslin(kappa,3)[[2]],choiceprobslin(kappa,4)[[2]]),ncol=4)

#find likelihood of observed data
like_wC=const_wC+sum(data_a1_wC*log(cprobvec_wC[1,]*(1-err)+err/3)+data_a2_wC*log(cprobvec_wC[3,]*(1-err)+err/3)+data_b1_wC*log(cprobvec_wC[2,]*(1-err)+err/3)+data_b2_wC*log(cprobvec_wC[4,]*(1-err)+err/3)+data_c1_wC*log((1-cprobvec_wC[1,]-cprobvec_wC[2,])*(1-err)+err/3)+data_c2_wC*log((1-cprobvec_wC[3,]-cprobvec_wC[4,])*(1-err)+err/3))
like_noC=const_noC+sum(data_a1_noC*log(cprobvec_noC[1,]*(1-err)+err/2)+data_a2_noC*log(cprobvec_noC[2,]*(1-err)+err/2)+data_b1_noC*log((1-cprobvec_noC[1,])*(1-err)+err/2)+data_b2_noC*log((1-cprobvec_noC[2,])*(1-err)+err/2))

#combine likes
return(as.numeric(-like_wC-like_noC))
}

if(indii %in% c(10,12,19)){
modellin=mle(likelihoodlin,start=list(kappa=10000,err=0.5), fixed=list(kappa=1000), method="L-BFGS-B",lower=c(0.001,0.001),upper=c(10000,0.999))
}else{
modellin=mle(likelihoodlin,start=list(kappa=7,err=0.1), method="L-BFGS-B",lower=c(0.001,0.001),upper=c(10000,0.999))
}


kappalin=coef(modellin)[1]
errlin=coef(modellin)[2]

lllin=logLik(modellin)[1]

#############################################
#Power Mutual Information

choiceprobspow=function(kappa,rho,treatment){
b1pay=b1vec[treatment]
b2pay=b2vec[treatment]
#Without C
#Information cost function
#cprobs in format a1,a2
cost_noC=function(cprobs){
overallcost=1/2*(cprobs[1]*log(cprobs[1]/(cprobs[1]+cprobs[2]))+(1-cprobs[1])*log((1-cprobs[1])/(2-cprobs[1]-cprobs[2]))+cprobs[2]*log(cprobs[2]/(cprobs[1]+cprobs[2]))+(1-cprobs[2])*log((1-cprobs[2])/(2-cprobs[1]-cprobs[2])))-2*0.5*log(0.5)
return(kappa*overallcost^rho)
}

#Optimization target
negpayoff_noC=function(cprobs){
return(-1*(r*(50/2*(cprobs[1]+cprobs[2])+b1pay/2*(1-cprobs[1])+b2pay/2*(1-cprobs[2]))-cost_noC(cprobs)))
}

#Find optimal accuracies
optimout_noC=constrOptim(rep(0.5,2),negpayoff_noC, grad=NULL,ui=conmat_noC,ci=convec_noC)
cprobvec_noC=optimout_noC$par

#With C
#Information cost function
#cprobs in format a1,b1,a2,b2
cost_wC=function(cprobs){
overallcost=1/2*(cprobs[1]*log(cprobs[1]/(cprobs[1]+cprobs[3]))+cprobs[2]*log(cprobs[2]/(cprobs[2]+cprobs[4]))+(1-cprobs[1]-cprobs[2])*log((1-cprobs[1]-cprobs[2])/(2-sum(cprobs)))+cprobs[3]*log(cprobs[3]/(cprobs[1]+cprobs[3]))+cprobs[4]*log(cprobs[4]/(cprobs[2]+cprobs[4]))+(1-cprobs[3]-cprobs[4])*log((1-cprobs[3]-cprobs[4])/(2-sum(cprobs))))-2*0.5*log(0.5)
return(kappa*overallcost^rho)
}

#Optimization target
negpayoff_wC=function(cprobs){
return(-1*(r*(50/2*(cprobs[1]+cprobs[3])+100/2*(1-cprobs[1]-cprobs[2])+b1pay/2*(cprobs[2])+b2pay/2*(cprobs[4]))-cost_wC(cprobs)))
}

#Find optimal accuracies
optimout_wC=constrOptim(rep(0.2,4),negpayoff_wC, grad=NULL,ui=conmat_wC, ci=convec_wC,)
cprobvec_wC=optimout_wC$par

return(list(cprobvec_noC,cprobvec_wC))

}

likelihoodpow=function(kappa,rho,err){
cprobvec_noC=matrix(c(choiceprobspow(kappa,rho,1)[[1]],choiceprobspow(kappa,rho,2)[[1]],choiceprobspow(kappa,rho,3)[[1]],choiceprobspow(kappa,rho,4)[[1]]),ncol=4)
cprobvec_wC=matrix(c(choiceprobspow(kappa,rho,1)[[2]],choiceprobspow(kappa,rho,2)[[2]],choiceprobspow(kappa,rho,3)[[2]],choiceprobspow(kappa,rho,4)[[2]]),ncol=4)

#find likelihood of observed data
like_wC=const_wC+sum(data_a1_wC*log(cprobvec_wC[1,]*(1-err)+err/3)+data_a2_wC*log(cprobvec_wC[3,]*(1-err)+err/3)+data_b1_wC*log(cprobvec_wC[2,]*(1-err)+err/3)+data_b2_wC*log(cprobvec_wC[4,]*(1-err)+err/3)+data_c1_wC*log((1-cprobvec_wC[1,]-cprobvec_wC[2,])*(1-err)+err/3)+data_c2_wC*log((1-cprobvec_wC[3,]-cprobvec_wC[4,])*(1-err)+err/3))
like_noC=const_noC+sum(data_a1_noC*log(cprobvec_noC[1,]*(1-err)+err/2)+data_a2_noC*log(cprobvec_noC[2,]*(1-err)+err/2)+data_b1_noC*log((1-cprobvec_noC[1,])*(1-err)+err/2)+data_b2_noC*log((1-cprobvec_noC[2,])*(1-err)+err/2))

#combine likes
return(-like_wC-like_noC)
}

if(indii %in% c(1,2,10,12,16,28)){
modelpow=mle(likelihoodpow,start=list(kappa=100,rho=1,err=0.1), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}else{
if(indii %in% c(13)){
modelpow=mle(likelihoodpow,start=list(kappa=1,rho=5,err=0.5), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}else{
if(indii %in% c(17,19,24)){
modelpow=mle(likelihoodpow,start=list(kappa=100,rho=1,err=0.9), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}else{
modelpow=mle(likelihoodpow,start=list(kappa=10,rho=1,err=0.4), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}
}
}


kappapow=coef(modelpow)[1]
rhopow=coef(modelpow)[2]
errpow=coef(modelpow)[3]

llpow=logLik(modelpow)[1]


#####################################################
#General Entropy

choiceprobsgen=function(kappa,rho,treatment){
b1pay=b1vec[treatment]
b2pay=b2vec[treatment]
#Without C
#Information cost function
#cprobs in format a1,a2
cost_noC=function(cprobs){
if(rho==1){
overallcost=kappa/2*(cprobs[1]*log(cprobs[1]/(cprobs[1]+cprobs[2]))+(1-cprobs[1])*log((1-cprobs[1])/(2-cprobs[1]-cprobs[2]))+cprobs[2]*log(cprobs[2]/(cprobs[1]+cprobs[2]))+(1-cprobs[2])*log((1-cprobs[2])/(2-cprobs[1]-cprobs[2])))
}
if(rho==2){
overallcost=kappa*(-1/2*1/2*(log(cprobs[1]/(cprobs[1]+cprobs[2]))+log((1-cprobs[1])/(2-cprobs[1]-cprobs[2]))+log(cprobs[2]/(cprobs[1]+cprobs[2]))+log((1-cprobs[2])/(2-cprobs[1]-cprobs[2])))-2*log(2))
}
if(rho!=1&rho!=2){
overallcost=kappa*(2^(1-rho)/((rho-2)*(rho-1))*1/2*((cprobs[1]/(cprobs[1]+cprobs[2]))^(2-rho)+(cprobs[2]/(cprobs[1]+cprobs[2]))^(2-rho)+((1-cprobs[1])/(2-cprobs[1]-cprobs[2]))^(2-rho)+((1-cprobs[2])/(2-cprobs[1]-cprobs[2]))^(2-rho))-1/((rho-2)*(rho-1))-log(2))
}
return(overallcost)
}

#Optimization target
negpayoff_noC=function(cprobs){
return(-1*(r*(50/2*(cprobs[1]+cprobs[2])+b1pay/2*(1-cprobs[1])+b2pay/2*(1-cprobs[2]))-cost_noC(cprobs)))
}

#Find optimal accuracies
optimout_noC=constrOptim(rep(0.5,2),negpayoff_noC, grad=NULL,ui=conmat_noC,ci=convec_noC)
cprobvec_noC=optimout_noC$par

#With C
#Information cost function
#cprobs in format a1,b1,a2,b2
cost_wC=function(cprobs){
if(rho==1){
overallcost=kappa/2*(cprobs[1]*log(cprobs[1]/(cprobs[1]+cprobs[3]))+cprobs[2]*log(cprobs[2]/(cprobs[2]+cprobs[4]))+(1-cprobs[1]-cprobs[2])*log((1-cprobs[1]-cprobs[2])/(2-sum(cprobs)))+cprobs[3]*log(cprobs[3]/(cprobs[1]+cprobs[3]))+cprobs[4]*log(cprobs[4]/(cprobs[2]+cprobs[4]))+(1-cprobs[3]-cprobs[4])*log((1-cprobs[3]-cprobs[4])/(2-sum(cprobs))))
}
if(rho==2){
overallcost=kappa/2*(log(cprobs[1]/(cprobs[1]+cprobs[3]))+log(cprobs[2]/(cprobs[2]+cprobs[4]))+log((1-cprobs[1]-cprobs[2])/(2-sum(cprobs)))+log(cprobs[3]/(cprobs[1]+cprobs[3]))+log(cprobs[4]/(cprobs[2]+cprobs[4]))+log((1-cprobs[3]-cprobs[4])/(2-sum(cprobs))))
}
if(rho!=1&rho!=2){
overallcost=kappa/(2*(rho-2)*(rho-1))*((cprobs[1]/(cprobs[1]+cprobs[3]))^(2-rho)+(cprobs[3]/(cprobs[1]+cprobs[3]))^(2-rho)-1+(cprobs[2]/(cprobs[2]+cprobs[4]))^(2-rho)+(cprobs[4]/(cprobs[2]+cprobs[4]))^(2-rho)-1+((1-cprobs[1]-cprobs[2])/(2-sum(cprobs)))^(2-rho)+((1-cprobs[3]-cprobs[4])/(2-sum(cprobs)))^(2-rho)-1 )
}
return(overallcost)
}

#Optimization target
negpayoff_wC=function(cprobs){
return(-1*(r*(50/2*(cprobs[1]+cprobs[3])+100/2*(1-cprobs[1]-cprobs[2])+b1pay/2*(cprobs[2])+b2pay/2*(cprobs[4]))-cost_wC(cprobs)))
}

#Find optimal accuracies
optimout_wC=constrOptim(rep(0.3,4),negpayoff_wC, grad=NULL,ui=conmat_wC, ci=convec_wC,)
cprobvec_wC=optimout_wC$par

return(list(cprobvec_noC,cprobvec_wC))

}

likelihoodgen=function(kappa,rho,err){
cprobvec_noC=matrix(c(choiceprobsgen(kappa,rho,1)[[1]],choiceprobsgen(kappa,rho,2)[[1]],choiceprobsgen(kappa,rho,3)[[1]],choiceprobsgen(kappa,rho,4)[[1]]),ncol=4)
cprobvec_wC=matrix(c(choiceprobsgen(kappa,rho,1)[[2]],choiceprobsgen(kappa,rho,2)[[2]],choiceprobsgen(kappa,rho,3)[[2]],choiceprobsgen(kappa,rho,4)[[2]]),ncol=4)

#find likelihood of observed data
like_wC=const_wC++sum(data_a1_wC*log(cprobvec_wC[1,]*(1-err)+err/3)+data_a2_wC*log(cprobvec_wC[3,]*(1-err)+err/3)+data_b1_wC*log(cprobvec_wC[2,]*(1-err)+err/3)+data_b2_wC*log(cprobvec_wC[4,]*(1-err)+err/3)+data_c1_wC*log((1-cprobvec_wC[1,]-cprobvec_wC[2,])*(1-err)+err/3)+data_c2_wC*log((1-cprobvec_wC[3,]-cprobvec_wC[4,])*(1-err)+err/3))
like_noC=const_noC+sum(data_a1_noC*log(cprobvec_noC[1,]*(1-err)+err/2)+data_a2_noC*log(cprobvec_noC[2,]*(1-err)+err/2)+data_b1_noC*log((1-cprobvec_noC[1,])*(1-err)+err/2)+data_b2_noC*log((1-cprobvec_noC[2,])*(1-err)+err/2))

#combine likes
return(-like_wC-like_noC)
}

if(indii %in% c(16)){
modelgen=mle(likelihoodgen,start=list(kappa=7.7,rho=7.9,err=0.1), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}else{
if(indii %in% c(20)){
modelgen=mle(likelihoodgen,start=list(kappa=4,rho=0.01,err=0.3), fixed=list(rho=0.01), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}else{
if(indii %in% c(15,17,19)){
modelgen=mle(likelihoodgen,start=list(kappa=7,rho=1,err=0.01), fixed=list(rho=1), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}else{
modelgen=mle(likelihoodgen,start=list(kappa=0.7,rho=0.9,err=0.1), method="L-BFGS-B",lower=c(0.001,0.01,0.001),upper=c(10000,100,0.999))
}
}
}

kappagen=coef(modelgen)[1]
rhogen=coef(modelgen)[2]
errgen=coef(modelgen)[3]
llgen=logLik(modelgen)[1]

##################################
#Inattentive Model

likelihoodinat=function(err){
cprobvec_wC=matrix(rep(c(0.5,0,0.5,0),4),ncol=4)
cprobvec_noC=matrix(rep(c(1,1),4),ncol=4)

#find likelihood of observed data
like_wC=const_wC++sum(data_a1_wC*log(cprobvec_wC[1,]*(1-err)+err/3)+data_a2_wC*log(cprobvec_wC[3,]*(1-err)+err/3)+data_b1_wC*log(cprobvec_wC[2,]*(1-err)+err/3)+data_b2_wC*log(cprobvec_wC[4,]*(1-err)+err/3)+data_c1_wC*log((1-cprobvec_wC[1,]-cprobvec_wC[2,])*(1-err)+err/3)+data_c2_wC*log((1-cprobvec_wC[3,]-cprobvec_wC[4,])*(1-err)+err/3))
like_noC=const_noC+sum(data_a1_noC*log(cprobvec_noC[1,]*(1-err)+err/2)+data_a2_noC*log(cprobvec_noC[2,]*(1-err)+err/2)+data_b1_noC*log((1-cprobvec_noC[1,])*(1-err)+err/2)+data_b2_noC*log((1-cprobvec_noC[2,])*(1-err)+err/2))

#combine likes
return(-like_wC-like_noC)
}

if(indii %in% c(1,19,22,24)){
modelinat=mle(likelihoodinat,start=list(err=0.01),  method="Brent",lower=c(0.01),upper=c(0.99))
}else{
modelinat=mle(likelihoodinat,start=list(err=0.001),  method="Brent",lower=c(0.0001),upper=c(0.99))
}

errinat=coef(modelinat)[1]
llinat=logLik(modelinat)[1]

addframe=data.frame(uid=uidlist[indii],likelin=lllin,errlin=errlin,likepow=llpow,errpow=errpow,likegen=llgen,errgen=errgen,likeinat=llinat,errinat=errinat)

indioutframe=rbind(indioutframe,addframe)

print(indii)
}

indioutframe$AIClin=indioutframe$likelin*-2+4
indioutframe$AICpow=indioutframe$likepow*-2+6
indioutframe$AICgen=indioutframe$likegen*-2+6
indioutframe$AICinat=indioutframe$likeinat*-2+2


indioutframe$minAIC=0
for(i in 1:length(indioutframe$AIClin)){
indioutframe$minAIC[i]=min(c(indioutframe$AIClin[i],indioutframe$AICpow[i],indioutframe$AICgen[i],indioutframe$AICinat[i]))
}

indioutframe$minAIClin=(indioutframe$AIClin==indioutframe$minAIC)
indioutframe$minAICpow=(indioutframe$AICpow==indioutframe$minAIC)
indioutframe$minAICgen=(indioutframe$AICgen==indioutframe$minAIC)
indioutframe$minAICinat=(indioutframe$AICinat==indioutframe$minAIC)

sum(indioutframe$minAIClin)/length(indioutframe$AIClin)
sum(indioutframe$minAICpow)/length(indioutframe$AIClin)
sum(indioutframe$minAICgen)/length(indioutframe$AIClin)
sum(indioutframe$minAICinat)/length(indioutframe$AIClin)



write.table(indioutframe, "C:/Users/neligh/Dropbox (Chapman)/Dean Stuff/Full Data Reports/experiment_1_individuals.csv", sep = ",", col.names = T, append = F)



